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1.  INTRODUCTION 


3.  ADAPTIVE  MULTIGRID  METHOD 


Accurate  prediction  of  hurricane  tracks  may  require  re¬ 
solving  the  flow  both  within  and  around  the  storm.  Since 
the  spatial  scales  in  these  two  regions  differ  substantially, 
uniform  resolution  is  inherently  inefficient:  the  grid  should 
be  refined  only  near  the  storm.  This  paper  describes  the 
performance  of  a  barotropic  model  with  an  adaptive  multi¬ 
grid  scheme  which  automatically  refines  the  mesh  around 
the  storm  as  it  moves. 

2.  MODEL  DESCRIPTION 

We  formulate  the  model  on  a  section  of  the  sphere 
using  a  Mercator  projection  (true  at  latitude  (j>  —  (pc)-  The 
model  consists  of  the  modified  barotropic  vorticity  equation 

^+rn^J{ip,0+Pm^  (1) 

where  the  relative  vorticity  C  and  streamfunction  are  re¬ 
lated  by 


Here  jdx^  'S  the  Jacobian 

of  (V?,C)  with  respect  to  (a;,y),  P  —  2Q  cos  0/a  (with 
a  and  Q  the  radius  and  rotation  rate  of  the  earth),  and 
m  ~  cos  (pc/ cos(p  is  the  map  factor.  There  are  two  quasi¬ 
physical  parameters:  the  diffusion  coefficient  z/,  and  the 
parameter  7  (inverse  of  the  effective  Rossby  radius)  which 
helps  prevent  retrogression  of  ultralong  Rossby  waves.  The 
model  domain  is  a  rectangle  in  x  and  y  centered  at  (x,  y)  = 
(0, 0),  where  (A,  p)  —  (Ac,  pc)-  At  the  boundaries  we  spec¬ 
ify  the  streamfunction  p  (and  thus  the  normal  component 
of ithe ’Velocity);  where  there  is  inflow,  we  also  specify  the 
vorticity 

]'•  The  space  discretization  uses  second-order  finite  dif¬ 
ferences  on  uniform  rectangular  grids  (as  explained  below), 
afiproxi mating  the  advection  terms  by  the  Arakawa  Jaco¬ 
bian.  The  time  discretization  uses  the  classical  fourth-order 
Runge-Kutta  (RK4)  scheme;  this  is  highly  accurate,  allows 
relatively  large  time  steps  for  stability,  and — since  it  is  a 
o'ljle-step  scheme — is  easy  to  implement  and  has  no  spurious 
computational  modes.  These  discretizations  are  embedded 
in  an  adaptive  method  which  includes  local  refinements  in 
both, space  and  time  as  detailed  in  the  next  section. 


*  Corresponding  author  address:  Scott  R.  Fulton,  Dept,  of 
Mathematics  and  Computer  Science,  Clarkson  University, 
Potsdam,  NY  13699-5815;  email:  fulton@clarkson.edu. 


To  adapt  the  resolution  near  the  storm  we  superimpose 
nested  uniform  grids  with  different  mesh  sizes.  The  base 
grid  Gi  covers  the  the  entire  computational  domain,  while 
successively  finer  patches  G2,  Gs,  ...  are  strictly  nested 
as  shown  in  Fig.  1.  The  basic  time  stepping  algorithm  is 
similar  to  that  used  in  most  nested  models  (e.g.,  Berger 
and  Oliger,  1984).  Specifically,  with  two  grids  (coarse  and 
fine)  one  full  time  step  is  executed  as  follows: 

1.  One  step  (length  At)  on  the  coarse  grid, 

2.  Two  steps  (length  At/2)  on  the  fine  grid,  with  bound¬ 
ary  values  interpolated  from  the  coarse  grid  in  space 
and  time, 

3.  Transfer  the  fine-grid  solution  to  the  coarse  grid  where 
they  overlap. 

This  algorithm  generalizes  recursively  to  more  than  two 
grids.  Each  time  step  consists  of  four  RK4  stages,  each 
of  which  involves  predicting  ^  from  (1)  and  diagnosing  p 
from  (2). 

To  solve  for  p  on  any  computational  grid  Gi  we  use  a 
multigrid  method.  With  local  time  stepping  as  above,  the 
next  coarser  grid  Gi-i  may  be  at  a  different  time  t,  so  the 
multigrid  method  usually  requires  using  an  additional  set  of 
coarse  grids  covering  the  same  domain  as  If  these  “local 
coarse  grids"  are  used  at  all  times,  the  resulting  method 
includes  only  one-way  interaction  between  the  grids:  the 
fine  grid  uses  boundary  values  from  the  coarse  grid,  but  Its 
effect  is  only  felt  on  the  coarse  grid  through  the  solution 
transfer  in  the  overlap  region  at  the  end  of  the  time  step. 
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Figure  1.  Computational  grids  {h  =  128,  64,  32,  16  km). 
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ACCURACY  VS.  EFFICIENCY 


The  local  coarse  grids  must  be  used  in  the  first  time 
step  on  the  fine  grid  Gi  in  the  Berger-Oliger  algorithm. 
However,  in  the  final  RK4  stage  of  the  second  fine-grid  time 
step — when  the  next  coarser  computational  grid  Gi-i  is  at 
the  same  time  level — the  full  computational  grid  Gi-i  may 
be  used  instead  of  the  first  local  coarse  grid  in  the  multigrid 
cycling.  In  this  way,  fine-grid  information  can  affect  the  en¬ 
tire  coarse  grid  immediately  through  the  relaxation  process, 
so  the  multigrid  algorithm  automatically  achieves  two-way 
interaction  between  the  computational  grids.  Note  that 
this  approach  requires  Full  Approximation  Scheme  (FAS) 
processing,  since  each  grid  patch  covers  a  smaller  domain 
than  the  previous  one.  Specifically,  if  we  denote  the  dis¬ 
crete  approximation  to  (2)  on  the  fine  grid  Gi  with  mesh 
size  h  =  hi 

=  F^  (3) 

then  the  corresponding  equation  to  be  solved  in  the  overlap 
region  on  the  coarse  grid  Gi-\  with  mesh  size  2h  =  /i/_i 
is 
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Figure  2.  Accuracy  vs.  efficiency  (from  Fulton,  1997) 


lih^2h  ^  p2h  (4) 

and  the  corresponding  correction  to  the  fine-grid  solution 
is  given  by 

(5) 

Here,  represents  the  approximate  solution  for  xp  on 
the  fine  grid,  the  fine-to-coarse  transfer  operators  and 
represent  injection  and  full  weighting,  respectively,  and 
the  coarse- to- fine  transfer  operator  I^h  represents  bilinear 
interpolation.  In  the  region  of  the  coarse  grid  not  cov¬ 
ered  by  the  fine  grid,  the  FAS  equation  (4)  is  replaced  by 
_  p2h  correction  (5)  is  needed. 

4.  SELF-ADAPTIVITY 

In  previous  versions  of  the  model  (e.g.,  Fulton,  1997) 
the  grid  sizes  were  fixed  in  advance,  and  simply  moved  when 
necessary  to  keep  the  location  of  the  vorticity  maximum 
roughly  centered  on  each  patch.  Figure  2  shows  sample 
results  quantifying  the  gains  in  efficiency  available  with  the 
mesh  refinement  strategy  given  above.  Each  point  gives  the 
mean  track  error  over  72  h  for  a  single  model  run;  the  label 
gives  the  finest  mesh  size  (in  km)  and  size  of  the  patches,  if 
any  (side  length  as  a  fraction  of  the  full  domain  width,  with 
A=l/2,  B=3/8,  C=l/4,  D=3/16,  and  E=l/8).  For  exam¬ 
ple,  the  grids  in  Fig.  1  are  for  the  run  128BCD.  The  gain  in 
efficiency  due  to  mesh  refinement  varies  significantly,  with 
the  best  combinations  of  patches  giving  speedups  in  the 
range  10-20.  Since  the  optimal  combination  of  patch  sizes 
will  depend  on  both  the  vortex  and  the  surrounding  flow 
and  may  vary  as  the  solution  evolves,  a  fully  self-adaptive 
procedure  would  be  advantageous.  Here  we  introduce  such 
a  method,  using  truncation  error  estimates  to  decide  auto¬ 
matically  where  to  refine  or  coarsen  the  mesh. 

The  necessary  truncation  error  may  be  estimated  es¬ 
sentially  for  free  using  FAS  processing  as  described  above. 
Specifically,  the  relative  truncation  error 

rl'*  =  F^'*  -  f'*  (6) 


provides  a  fourth-order  approximation  to  the  difference  in 
truncation  error  between  grids  h  and  2h.  Thus,  the  trun¬ 
cation  error  on  grid  h  can  be  approximated  with  accuracy 
0{h^)  as  "^^is  estimate  is  computed  at  each 

point  on  grid  h  at  the  end  of  the  time  step.  Where  is 
small,  the  current  fine  grid  h  is  no  longer  needed,  so  its  ex¬ 
tent  may  be  reduced.  Conversely,  where  is  large,  further 
refinement  may  be  beneficial,  so  these  points  are  flagged 
and  a  new  finer  grid  h/2  enclosing  them  is  introduced.  The 
determination  of  what  is  "smaH"  or  “large"  follows  the  al¬ 
gorithm  outlined  in  Brandt  (1977,  sec.  8),  which  uses  a 
parameter  A  to  control  the  trade-off  between  increased  ac¬ 
curacy  and  increased  work.  Note  that  since  and 
are  already  computed  in  the  solution  process,  computing 
the  truncation  error  estimate  requires  very  little  additional 
work.  Also,  the  estimated  truncation  error  can  be  used  in 
an  extrapolation  method  to  obtain  fourth-order  accuracy 
from  the  second-order  discretization.  We  will  present  re¬ 
sults  at  the  conference  quantifying  the  performance  of  the 
self-adaptive  method  outlined  above. 
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